function distance2=distance2(x)
[ca_data] = xlsread('CA_mexico_data');
ca_data=ca_data./100;
SS_data=zeros(1,length(ca_data));
SS_data(ca_data>(mean(ca_data)+2*std(ca_data)) )=1;
Proba_data=mean(SS_data);
auto_data=ar(ca_data,1);
ar_data=[auto_data.A(2),auto_data.NoiseVariance];

% % VECTOR OF DATA PARAMETERS
global data_par Tradables_share Init_Par

%data_par=[-0.319567456685736,0.055,0.02257994] ; %Debt; Probability of Crisis; Standard deviation of CA/Y
data_par=[-0.319567456685736,Proba_data,ar_data(2),ar_data(1)] ; %Debt; Probability of Crisis; Standard deviation of CA/Y
Tradables_share=0.447147283578164; %Share of Production of tradable goods to GDP

Weight=[1/abs(data_par(1)),0,0,0;0,1/abs(data_par(2)),0,0;0,0,1/abs(data_par(3)),0;0,0,0,1/abs(data_par(4))];

%Vector of Initial Parameters
Par=[x(1),x(2),x(3),x(4)];


%%%%%First Set of Parameters%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%  Technical parameters  %%%%%%%%%%%%%%%%%%%%%%%%%
uptd=0.2; %Change to 0.1 if policy function fluctuate
outfreq=20; %Display frequency
iter_tol=600; % Iteration tolerance
tol=1e-5; %Numerical tolerance
 %Number of States
NB=100; %number of nodes.
bmin= -1;  %-1.11;
bmax=-0.;
NT      = 5;                % Number of grid points for tradables
NR      = 1;                % Number of grid points for nontradables
%NKAPPA  = 2;
NKAPPA  = 3;
NSS=NR*NT*NKAPPA;

%Moments from the Data
Debt_Data=data_par(1);
ProbabilityOFSS=data_par(2);
Std_CA=data_par(3);
ShareofT=Tradables_share;

%CALIBRATED PARAMETERS
beta=Par(1);
meank=Par(2);
sigmak=Par(3);

rhok=Par(4);

%PARAMETERS NOT CALIBRATED
sigma=2;
ita=(1/0.83)-1;
rbar=0.04;
yN=1;

%AR(1) process for income
rho=0.460247378300448 ;
sigmay=0.032352010489055; 
 
[Z,Zprob] = tauchenhussey(NT,0,rho,sigmay,sigmay); % Row are today, column are tomorrow
[EiVec, EiVal]=eig(Zprob);
EiVec1=EiVec(:,1);
num1=sum(EiVec1);
EiVec1=EiVec1/num1; % The probability of Each of the 5 States of YT, in  this case constant




%AR(1) process for Kappas
%rho_mod=rho+change;
%rhok=rho_mod;
%rhok=rhok0

[NewK,NewKprob]= tauchenhussey(NKAPPA,meank,rhok,sigmak,sigmak);


%Prob_kappa=[ 0.8 0.2; 0.2 0.8];

%Prob=kron(Prob_kappa,Zprob);   % check this is the right
Prob=kron(NewKprob,Zprob);   

%kap_values= [0.32 0.5]';

yT      = exp(Z);
Sr     =  ones(NB,NR)*exp(rbar);
 

YT=repmat(yT,NKAPPA,NB )';
SR=repmat(Sr,1,NSS );
%KAPPAS=kron(kap_values,ones(NT,NB))';
KAPPAS=kron(NewK,ones(NT,NB))';

%Calibration of OMEGA
%omega=0.3005; %Calibrated Outside of The Model
RatioOmega=mean((yT.*(1-ShareofT)).*(ShareofT.*(yT-Debt_Data*(exp(rbar)-1)).^(1+ita)).^(-1));
omega=1/(1+RatioOmega);



%%%%%%%%%%%%%%%%%%%%%%%%%%%Second Set of Parameters%%%%%%%%%%%%%%%%%%%%%%%%


B=zeros(NB,1);
B(1)=bmin;
 
for i=2:NB;
    B(i)=B(i-1)+(bmax-B(i-1))/((NB-i+1)^1.05); 
end
 
 


%Q=1./SR;
q=1/exp(rbar);

b=repmat(B,1,NSS);
yn=ones(NB,NSS);
 
bp=  b;
%c= b./SR  +YT  -bp;
c= b*q  +YT  -bp;
%price=real((1-omega)/omega*c.^(1+ita));
price=real((1-omega)/omega*c.^(1+ita));

c=ones(NB,NSS);
EE=zeros(NB,NSS);
cbind=ones(NB,NSS);
%bpmax=-kappa*(price+YT);
bpmax=-(1/q)*KAPPAS.*(price+YT);

% mup=ones(NB,NSS);
% emu=zeros(NB,NSS);
%mg_utilbind=zeros(NB,NSS);
c_p=ones(NB,NSS,NSS);
c_tom=ones(1,NSS);
totalcp=ones(1,NSS);
mu_p=ones(1,NSS);
emu=0;


cbind=c;
 
uptd=0.2; %Change to 0.1 if policy function fluctuate
outfreq=20; %Display frequency
iter_tol=600; % Iteration tolerance
tol=1e-5; %Numerical tolerance

%%%%%%%%%%%%%%%%%%DECENTRALIZED EQ%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear kappa
load A_de

iter=0;
d2=100;
options = optimoptions('fsolve','Display','off');
disp('DE Iter      Norm');

  tol=0.00001
   iter=1
   
while d2>tol && iter <iter_tol    
    
oldp=price;
oldc=c;
oldbp=bp;

c_p=interp1(B,c,bp,'linear','extrap'); %This is 3D Array with consumptions tomorrow in today's grid
%mg_utilbind=(omega*cbind.^(-ita)+(1-omega)*yn.^(-ita)).^(sigma/ita -1/ita -1 )*omega.*cbind.^(-ita-1);
%  totalc=(omega*c.^(-ita)+(1-omega)*yn.^(-ita)); 
%  mup=real(omega*totalc.^(sigma/ita).*(totalc.^(-1/ita-1)).*(c.^(-ita-1)));  %mup(b,y)
% 

%   parfor i=1:NB
%        for j=1:NSS
%             emu(i,j) =beta*exp(rbar)*interp1(B,mup,bp(i,j),'linear','extrap')*Prob(j,:)' %EMU is expected marginal utility tomorrow in today's grid. 
%        end
%   end



parfor i=1:NB
    for j=1:NSS
         c_tom=squeeze(c_p(i,j,:))'; % This is a vector with consumption tomorrow in each of the possible 15 states when today state is (i,j)
         totalcp=(omega*c_tom.^(-ita)+(1-omega)*yn(1,:).^(-ita));
         mu_p=real(omega*totalcp.^(sigma/ita).*(totalcp.^(-1/ita-1)).*(c_tom.^(-ita-1)));
         emu=beta*exp(rbar)*mu_p*Prob(j,:)'; %Expected Marginal Utility at state (i,j)
         EE(i,j)=(omega*cbind(i,j)^(-ita)+1-omega)^(sigma/ita -1/ita -1 )*omega*cbind(i,j)^(-ita-1)- emu;
%          EE(i,j)=(omega*cbind(i,j)^(-ita)+1-omega)^(sigma/ita -1/ita -1 )*omega*cbind(i,j)^(-ita-1)- emu(i,j);
        if EE(i,j)>1e-5
            bp(i,j)=bpmax(i,j);
            c(i,j)=cbind(i,j);
       
        else                            
            %f = @(cc) (omega*cc^(-ita)+1-omega)^(sigma/ita-1/ita-1)*omega*cc^(-ita-1)-emu(i,j);
            f = @(cc) (omega*cc^(-ita)+1-omega)^(sigma/ita-1/ita-1)*omega*cc^(-ita-1)-emu;
            x0=c(i,j);
            [c(i,j), EE(i,j)]=fzero(f,x0);           
            bp(i,j)=max((1/q)*(YT(i,j)+b(i,j)-c(i,j)),bmin);
            bp(i,j)=min(bp(i,j),bmax);
        end
           
 
    end
end
 
price=real((1-omega)/omega*c.^(1+ita));

c=b+YT-max(q*bp,-KAPPAS.*(price+YT));  %qb'>=-kappa*(yt+yn*pn)
price=real((1-omega)/omega*c.^(1+ita));
bp=(1/q)*(b+YT-c);

d2=max([max(max(abs(c-oldc))),max(max(abs(bp-oldbp)))]);

%=====================Updating rule.Must be slow, important!==============
bp=uptd*bp+(1-uptd)*oldbp;
c=uptd*c+(1-uptd)*oldc;
price=uptd*price+(1-uptd)*oldp;
%========================================================================

bpmax=-(1/q)*KAPPAS.*(price.*yn+YT);
bpmax(bpmax>bmax)=bmax;
bpmax(bpmax<bmin)=bmin;
cbind=b+YT-q*bpmax;

 
cbind(cbind<0)=nan;


iter=iter+1;


D2(i,:)=d2;

metric(i)=d2;
if mod(iter, outfreq) == 0;
    fprintf('%d          %1.7f \n',iter,d2);
end
end
save de
clear  KAPPAS beta meank sigmak NewK NewKprob bpmax Prob Par beta0 meank0 sigmak0 iter rhok0 rhok
save A_de
clear all
load de

%%%%%%%%%%%%%%%%%%%%%%%%%%Simulation%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ISOM_01_model
%%%%%%%%%%%%% Parameters %%%%%%%%%%%%%%%%
%clear all
T=500000;
cut=100;  % Initial condition dependence
%load model_Brazil_k
%addpath('auxiliar')
clear kappa

%% 1. Simulation

[S_index,]=markov(Prob,T+cut+1,1,1:length(Prob));

    ySIM=YT(1,S_index)';  
    RSIM=SR(1,S_index)';  
    KAPPASIM=KAPPAS(1,S_index)';

bpSIM=zeros(T+cut,1);
%WfSIM=zeros(T+cut,1);
bindSIM=zeros(T+cut,1);
cSIM=zeros(T+cut,1);
pSIM=zeros(T+cut,1);

%bpSP_SIM=zeros(T+cut,1);
%bindSP_SIM=zeros(T+cut,1);
%cSP_SIM=zeros(T+cut,1);
%pSP_SIM=zeros(T+cut,1);

%Tau_SIM=zeros(T+cut,1);
%Tregion_SIM=zeros(T+cut,1);
 

bpSIM(1)=bmax+(bmin-bmax)/2;
%bpSP_SIM(1)=bmax+(bmin-bmax)/2;

for i=2:T+cut
     % de sim
    bpSIM(i)=interp1(B,bp(:,S_index(i)),bpSIM(i-1),'linear');
   cSIM(i)=bpSIM(i-1)+ySIM(i)-(1/(RSIM(i)))*bpSIM(i); 
 if isnan(bpSIM(i))==1
     %pause
 end
   % WfSIM(i)=interp1(B,Wf(:,S_index(i)),bpSIM(i-1),'linear');

    
         % sp sim
    %bpSP_SIM(i)=interp1(B,bpSP(:,S_index(i)),bpSP_SIM(i-1),'linear');
    %cSP_SIM(i)=bpSP_SIM(i-1)+ySIM(i)-(1/(RSIM(i)))*bpSP_SIM(i); 
    %Tau_SIM(i)=interp1(B,tau_sp(:,S_index(i)),bpSP_SIM(i-1),'linear');
end

pSIM=(1-omega)/omega*cSIM.^(1+ita); 
%pSP_SIM=(1-omega)/omega*cSP_SIM.^(1+ita);

bplim_SIM=-(RSIM).*  KAPPASIM.*(pSIM+ySIM);
%bplimSP_SIM=-(RSIM).*  KAPPASIM.*(pSP_SIM+ySIM);

%%%%% IMPORTANT TO USE TOLERANCE NOT <= %%%%%%%%%%%%
bindSIM((bpSIM-bplim_SIM)<=1e-5)=1;
%bindSP_SIM((bpSP_SIM-bplimSP_SIM)<=1e-5)=1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

bpSIM(bindSIM==1)=bplim_SIM(bindSIM==1);
%bpSP_SIM(bindSP_SIM==1)=bplimSP_SIM(bindSP_SIM==1);

%Tregion_SIM(Tau_SIM>1)=1;

CA_SIM=bpSIM(2:end)-bpSIM(1:end-1);  % >0 capital outflow
%CASP_SIM=bpSP_SIM(2:end)-bpSP_SIM(1:end-1);

% Eliminate Initial Condition Dependency
%S_index=S_index(cut+1:T+cut)';

bpSIM=bpSIM(cut+1:T+cut);
%bindSIM=bindSIM(cut+1:T+cut);
%bplim_SIM=bplim_SIM(cut+1:T+cut);
%cSIM=cSIM(cut+1:T+cut);
pSIM=pSIM(cut+1:T+cut);
CA_SIM=CA_SIM(cut:T+cut-1);

%bpSP_SIM=bpSP_SIM(cut+1:T+cut);
%bindSP_SIM=bindSP_SIM(cut+1:T+cut);
%bplimSP_SIM=bplimSP_SIM(cut+1:T+cut);
%cSP_SIM=cSP_SIM(cut+1:T+cut);
%pSP_SIM=pSP_SIM(cut+1:T+cut);
%CASP_SIM=CASP_SIM(cut:T+cut-1);

%WfSIM=WfSIM(cut+1:T+cut);
%Tau_SIM=Tau_SIM(cut+1:T+cut);
%Tregion_SIM=Tregion_SIM(cut+1:T+cut);
ySIM=ySIM(cut+1:T+cut);
%KAPPASIM=KAPPASIM(cut+1:T+cut);
%RSIM=RSIM(cut+1:T+cut);

Y_SIM=pSIM+ySIM;   % Income
%YSP_SIM=pSP_SIM+ySIM;

CA_SIM=CA_SIM./Y_SIM;  % Current Account
%CASP_SIM=CASP_SIM./YSP_SIM;

%CAchg_SIM=CA_SIM(2:end)-CA_SIM(1:end-1);
%CAchg_SIM=[0;CAchg_SIM];

Lev_SIM=bpSIM./Y_SIM;  % Leverage
%LevSP_SIM=bpSP_SIM./YSP_SIM;

%RER_SIM=(omega^(1/(1+ita))+(1-omega)^(1/(1+ita)).*pSIM.^(ita/(1+ita))).^((-1-ita)/ita); % Real Exchange Rate
%RERSP_SIM=(omega^(1/(1+ita))+(1-omega)^(1/(1+ita)).*pSP_SIM.^(ita/(1+ita))).^((-1-ita)/ita);


%MOMENTS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear SS
SS=zeros(T,1);
SS(CA_SIM>(mean(CA_SIM)+2*std(CA_SIM)) )=1; % Crisis is defined as current account goes 2 sd away and collateral constraint binds in decentralized economy
MM1=mean(Lev_SIM);
MM2=mean(SS);
auto_model=ar(CA_SIM,1);
MM4=auto_model.A(2);
MM3=auto_model.NoiseVariance;
moments=[MM1,MM2,MM3,MM4];

distance2 =(data_par-moments)*Weight*(data_par-moments)';
end